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A simulation and digital synthesis method of an oscillating 



The invention refers to a digital simulation method of a 
non-linear interaction between an excitation source and a wave 
in a resonator and may be applied, in particular, to the digital 
synthesis, in real time, of an oscillating phenomenon such as 
the sound emitted by a musical instrument operating more 
particularly with sustained oscillations, such as wind or rubbed 
string instrument. 

The phenomena of wave propagation and formation of the 
emitted sounds, in particular, by a musical instrument have 
been studied scientifically for a very long time. 

In particular, it is admitted, generally speaking, that a 
musical instrument includes, at least, one exciter, characterised 
by a non-linear characteristic, coupled possibly with certain 
linear elements (the reed, the lips, the bow, the hammer, etc..) 
and resonator elements, generally linear, where there is wave 
propagation as well as, generally, localised elements (for 
instance lateral bores or simple elements of the mass or spring 
type), generally linear as well. 

Similarly, a digital instrument capable of synthesising the 
sounds emitted by a musical instrument, is composed generally 
of three main elements, respectively a first element, to sense 
the gests of a musician and to transform them into 
signals/control parameters, a second element computing the 
signal in real time, a third element converting this series of 
numbers calculated into a sound signal by means of 
digital/analogue converters, amplifier, loudspeakers. 

The present invention concerns mainly the second real 
time calculation element of the signal. 

It is known that the digital simulation of a sound or, more 
generally, of an oscillating phenomenon, may be conducted by 
discretisation in the time domain of equations forming the 
mathematic representation of the physical phenomenon to be 
simulated. Such a model is always expressed in the form of a 



phenomenon 



system of equations with coupled partial derivations, linear or 
non-linear. 

The simulation then consists, generally speaking, in 
computing as quickly as possible the solution of the 
5 acoustic/mechanical model describing the operation of the 
instrument or, at least, approximations preserving its most 
important characteristics. 

Numerous methods exist to that effect and it is possible, 
in particular, to mention modal methods (which describe the 

10 resonator as a resonant filter comprised of a sum of elementary 
resonances), particular methods (which describe the medium 
wherein there is a wave propagation in the form of chains of the 
type mass-springs-dampers), or the digital methods for solving 
equations with partial derivations. 

15 However, real time sound synthesis is difficult to realise 

and consequently, for some years, other methods have been 
developed, based on a "signal processing" formalism of the 
propagation in both directions of the resonator of the 
instrument. One may quote, for instance the methods called 

20 "digital wave guide" or "digital wave filter". 

Generally speaking, to represent the propagation of a 
wave, one may use, in the simplest formalism, the well known 
d'Alembert's equation, which applies to longitudinal waves 
(acoustic for instance) as well as transversal waves (vibration of 

25 a string for instance). In particular, in the case of the 
propagation of an acoustic wave, the acoustic pressure in all 
points of the resonator of a wind instrument may be split into a 
sum of two waves of acoustic pressure, one propagating from 
the player to the horn, and the other from the horn to the player, 

30 which are called away-wave and return-wave. 

In practice, such propagation is expressed by a 
convolution equation (linear filtering), which yields the away 
wave (or return wave) at one point of the resonator at each time 
in relation to the away wave (or return wave) at another point at 

35 each instant. In the so-called Green formulation, which may be 




implemented in digital form, d'Alembert's equation specifies for 
instance that this linear filter, called Green core, is a pure 
delay, depending on the speed of propagation in the medium 
and on its length. 
5 In a synthesis model, these waves, respectively away and 

return waves, are represented by two signals corresponding 
respectively to both propagative solutions of the differential 
equation. 

Such a synthesis method is implemented, for instance in 
10 the document US-A-5, 332,862 which describes a synthesiser 
comprising generally speaking: 

- one non-linear part, simulating the exciter, to which two 
control parameters of the sound to be simulated are applied. 
These parameters are, in this case, the pressure of the player's 

is breath and the pressure of his lips on the reed or the 
mouthpiece, 

- a linear part, simulating the resonator, which receives a 
signal noted q 0 , representative of the away wave, emitted by the 
non-linear part, and which emits a signal noted q it 

20 representative of the return wave towards said part, 

- a means for creating the sound from the signals derived 
from the linear part and the non-linear part, 

- a digital/analogue converter generating the synthesised 
sound. 

25 Obviously, there are other types of synthesisers but, until 

now, all the methods involving a modelization the physical 
phenomena within the instrument were based on the 
decomposition of the vibration inside the resonator in terms of 
away wave and return wave variables. 

30 Still, it has appeared that such methods exhibited several 

shortcomings. 

First of all, when the acoustic resonator is formed, for 
instance, of several parts of cylindrical tubes of different 
diameters, the section change causes the generation of a 
35 transmitted wave and of a reflected wave at each interface. The 



phenomenon has been taken into account for a long time, for 
instance, within the framework of the modelization of the vocal 
conduit. 

This type of modelization, which is identical, in its 
5 approach, to the conventional theory of the geometrical optics, 
is also employed, for instance, in seismic-reflection, in order to 
describe the propagation of elastic waves in a multilayer 
ground. 

It is known, indeed that it is interesting, in all the cases 

10 when one or several waves propagate, to characterise an 
interface by a diffusion matrix, since it is thus possible to 
access the reflexions and transmissions of the different waves 
directly. However, the behaviour of this localised element often 
becomes difficult to grasp and to calculate, insofar as the 

15 interface and continuity equations are always written initially 
with physical quantities, for instance by expressing, at the 
interface, the continuity of the pressure or of the flow, of the 
strength or of the speed. 

It is hence often more advantageous to use "impedance" 

20 or "admittance" matrices, which link the physical quantities 
directly, as described in an article of J. Kergomard, "Calculation 
of discontinuities in waveguides using mode-matching method: 
an alternative to the scattering matrix approach" J. d'acoustique 
4 pp 111-138, (1991)". 

25 On the other hand, when there are localised elements 

other than interfaces in the instrument to be simulated, the 
"waveguide" method must be complemented by a "wave filter" 
type method describing such localised elements (such as mass, 
spring, dampers) for a correct connection between the various 

30 sub-systems. 

Similarly, when the acoustic resonator is composed for 
instance of a conical pipe, the waves moving from the player to 
the horn (away waves) and from the horn to the player (return 
waves) are different, which requires, either different 

35 modelization, by means of two linear filters corresponding to the 




Green cores describing the propagation in each direction, or 
approximating the cone by a succession of short length 
cylinders of different diameters. 

Besides, the sound produced by a musical instrument is 
5 not derived solely from the propagation of a wave through a 
resonator, regardless of the complexity of its geometry, but 
results from the non-linear coupling between said resonator and 
an excitation source. Such non-linear coupling is expressed 
physically between the physical quantities representing a cause 
10 (pressure in the acoustic case, strength in the mechanical case) 
and an effect (flow in the acoustic case, speed in the 
mechanical case), called Kirchhoff variables. In the acoustic 
case, this so-called Euler-Bernoulli physical law, a simplified 
version of the Navier-Stokes equations used in fluid mechanics, 
15 specifies that the acoustic pressure, at the reed or the lips of a 
wind instrument, is proportional, within one additive constant, to 
the square of the acoustic flow. The formulation of the waves in 
the resonator in the form of away waves and return waves 
requires therefore changing the variables enabling the non- 
20 linear coupling to be expressed, not in relation to the pressure- 
flow physical variables any longer, but in relation to these new 
away-wave and return-wave variables, as for instance, in the 
aforementioned document US-A-5, 332,862. Still, the variable 
change introduces an additional complexity in the synthesis 
25 method. It is thus that the use of iterative or tabulation methods 
has been recently suggested to calculate the solution of the 
non-linear system. 

Consequently, the synthesis methods used until now do 
not enable to simply express the non-linear coupling which 
30 exists between the excitation source and the resonator of the 
instrument and limit the physical parameterisation of the 
synthesis algorithms. 

The invention intends to remedy such shortcomings and to 
eliminate such limitations thanks to a new real-time simulation 
35 and synthesis method of an oscillating phenomenon, applicable 




especially, but without being limited thereto, to self-oscillating 
wind instruments. In particular, the invention refers to a 
simulation method enabling to take into account the physical 
process governing the operation of a real instrument and where 
5 the digital implementation may be particularly simple. 

Moreover, from a basic method applicable to instruments 
such as a clarinet, with a cylindrical resonator, the invention 
may be adapted to the simulation of other types of wind or 
string instruments. 

10 Besides, the invention is not limited to the simulation of 

musical instruments but may be applied, generally to real-time 
digital synthesis of all sorts of oscillating phenomena. 

Generally, the invention therefore relates to the simulation 
of a non-linear interaction between an excitation source and a 

15 wave in a resonator, by means of a digital calculation tool, from 
equations whose solution corresponds to the physical event of a 
phenomenon to be simulated. 

According to the invention, the phenomenon to be 
simulated is translated, at each time and at a given point of the 

20 resonator, by a linear relation between two variables 
representative of the effect and of the cause of said 
phenomenon, the impedance or admittance equation is 
transcribed directly in the form of a digital model enabling to 
realise a non-linear interaction between the two variables of the 

25 impedance or admittance relation. 

In this view, the model comprises, on the one hand, at 
least one linear part representing directly the so-called input 
impedance or admittance of the resonator, i.e. at the point 
where the non-linear interaction occurs and, on the other hand, 

30 one non-linear part modelization the role of the excitation 
source of the phenomenon to be simulated. 

In particular, for real-time digital synthesis of an 
oscillating phenomenon, the invention enables, from a system of 
equations between at least two variables representative of the 

35 behaviour of the resonator, to establish an expression of the 
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input impedance or admittance of the resonator in the form of a 
linear filter including delays, without any decomposition into 
away-return waves, in order to realise at least one linear part of 
the model which may be coupled with a non-linear loop involving 
5 the evolution of the non-linearity as expressed between the two 
variables of the impedance or admittance relation of the 
resonator. 

Particularly advantageously, such linear part of the model 
is composed of the sum of two elementary waveguides fulfilling 

10 a transfer function between the two variables of the impedance 
or admittance relation. 

According to another particularly advantageous 
characteristic, the model is driven by at least two parameters 
representative of the non-linear physical interaction between the 

15 source and the resonator, by means of a loop connecting the 
output to the input of the linear part and comprising a non-linear 
function playing the part of an excitation source for the 
resonator. 

Thus, contrary to the synthesis methods used 

20 conventionally, the method according to the invention does not 
involve away and return waves, but expresses directly and 
digitally the so-called impedance, linear relation between the 
cause and effect variables, i.e. pressure and flow in the 
acoustic case, strength and speed in the mechanical case. 

25 Naturally, such relation puts in evidence propagative 

elements, involving filters and delays, insofar as the physical 
phenomenon to be simulated is unchanged. 

Still, thanks to the method according to the invention, 
such relation is readily processable in digital form and may then 

30 be associated with the non-linear relation expressed physically 
between the same variables. 

As indicated above, the present invention concerns 
therefore essentially the modelization element of a digital 
instrument which, from parameters prepared by a control 

35 means, such as a gestural sensor operated by the player, 




computes in real time a signal liable to be transformed into a 
sound signal by a conversion element. 

In particular, for real-time synthesis, by physical 
modelization, of the sound of a musical instrument resulting 

5 from a non-linear coupling between the excitation source and 
the resonator, the invention enables to solve the system of 
equations representative of the phenomenon to be simulated by 
expressing directly and digitally the impedance or admittance 
linear relation between the cause and effect variables and by 

10 associating such linear relation in digital form with the non- 
linear relation between the same variables. Moreover, in the 
case of a resonator of a complex geometry, the former may be 
decomposed into successive elements, in order to combine the 
elementary linear relations corresponding respectively to each 

15 element of the resonator, in order to obtain an impedance or 
admittance corresponding to the geometry of the instrument. 

As already mentioned, the invention applies, in particular, 
to real-time synthesis of the sound produced by a wind 
instrument. In such a case, the two variables of the impedance 

20 relation are the acoustic pressure and flow at the input of the 
resonator. 

In the case of a cylindrical resonator having one open end, 
it is particularly advantageous to realise the linear part of the 
digital transcription model of the impedance equation in the 
25 form of a sum of two elementary waveguides having as an 
excitation source the flow at the input of the resonator and 
fulfilling the transfer function : 

7 (a) = Pe((o) = 1 exp(-2ik(co)L) 

eK } Ue(co) 1 + exp(-2ik(co)L) 1 + exp(-2ik(co)L) 



30 



wherein: 

- co is the wave angular frequency, 

- Ze(co) is the input impedance of the resonator, 




- Pe(co) and Ue(co) are the Fourier transforms of the 
dimensionless values of the pressure and of the flow at 
the input of the resonator, 

- k(co) is a function of the wave angular frequency which 
5 depends on the phenomenon to be simulated, 

- L is the length of the resonator. 

According to another characteristic, each of the two 
waveguides involves a filter having as a transfer function: 

10 - F(co) 2 = - exp(-2ik(co)L) 

and representing a two-way travel of a wave with the sign 
changing at the open end of the resonator, each waveguide 
corresponding to a term of the impedance equation. 

15 Such a model may advantageously be driven by the length 

of the resonator and at least two parameters representative of 
the non-linear physical interaction between the pressure and 
flow at the input of the resonator, by means of a loop 
connecting the output to the input of the linear part and 

20 comprising a non-linear function playing the part of an 
excitation source for the resonator. 

In particular, for real-time synthesis of the sounds to be 
simulated, a formulation is prepared, in the time domain, of the 
angular frequency response of the resonator, by approximation 

25 of the losses represented by the filter by means of an 
approximated digital filter. 

The invention covers other essential characteristics 
mentioned in the claims and referring, in particular, to the 
equations used by the digital signal calculation tool and which 

30 leads to waveguide models depending on the phenomenon to be 
simulated. 

Indeed, according to an essential characteristic of the 
invention, the method suggested for the simulation of a simple 
phenomenon such as the propagation of a wave in a cylindrical 
35 resonator, may be adapted in multiple ways for the simulation of 




more complex phenomena and, in particular, of diverse types of 
instruments. 

In the following description, we shall thus expose in detail 
the simulation method, the equations used and the model to be 
5 implemented for synthesising the sound of an instrument with a 
reed acoustic cylindrical resonator, such as a clarinet, and then 
certain adaptations for the simulation of other types of 
instruments. 

Figure 1 schematically represents a whole digital 
10 instrument for simulating a wind instrument, by the method 
according to the invention. 

Figure 2 provides two diagrams respectively representing, 
on the left, the transfer function, in Hertz, of a single-mode 
reed model and, on the right, the angular frequency response in 
is relation to the samples, with a sampling frequency of 44,100 
Hertz. 

Figure 3 is diagram per a calculation by combination of 
waveguides, representing the input impedance of a cylindrical 
resonator. 

20 Figure 4 provides two diagrams respectively representing, 

for a cylindrical resonator, at the top the input impedance in 
relation to the frequency expressed in Hertz and, at the bottom, 
the angular frequency response in relation to time, in seconds. 
Figure 5 is a calculation diagram of a simulation model of 

25 a cylindrical resonator reed-type instrument. 

Figure 6 provides two diagrams similar to Figure 3, 
representing respectively, for a resonator model computed 
according to the invention, at the top the approximated input 
impedance and at the bottom the approximated impulse 

30 response. 

Figure 7 provides two diagrams similar to Figure 2, 
respectively representing, for a reed model computed according 
to the invention, on the left the transfer function and on the right 
the angular frequency response. 




Figure 8a shows the variations, in relation to time 
expressed in seconds, of the internal acoustic pressure at the 
mouthpiece of a cylindrical resonator. 

Figures 8b and 8c are enlargements of attack and 
5 extinction transients. 

Figure 9 provides two diagrams respectively representing, 
on the left, the transfer function and, on the right, the angular 
frequency response, for a multimode reed model computed 
according to the invention. 
io Figure 10 provides two diagrams representing the 

spectrum of the external acoustic pressure, respectively, at the 
top, for a single-mode reed and, at the bottom, for a multiple 
mode reed. 

Figure 11 is a calculation diagram representing the 
15 impedance of a cylindrical resonator with terminal impedance. 

Figure 12 is a calculation diagram representing the 
impedance of a conical resonator. 

Figure 13 is a calculation diagram representing the 
impedance of a resonator for a wind instrument. 
20 Figure 14 is a general calculation diagram representing 

the impedance of a parallel combination of cylindrical 
resonators. 

Figure 15 provides two diagrams respectively 
representing, in the case of a string, at the top the exact 
25 admittance and at the bottom the approximated admittance, in 
relation to frequency expressed in Hertz. 

Figure 16 is a model of a digital instrument simulating a 
string instrument. 

Figure 17 shows, for a string struck, the time variations, at 
30 the top, of the speed" of the string at the contact point and, at 
the bottom, of the strength exerted by the hammer on the string. 

Figure 18 represents, for a string struck, the trajectory of 
the strength with time, in relation to relative displacement of the 
hammer with respect to the string. 



Figure 19 is a general simulation diagram of an instrument 
operating by non-linear coupling between an excitation source 
and a resonator. 

The invention will first of all be described in its application 

5 to a clarinet-type wind instrument. 

Figure 1 schematically represents a whole digital 
instrument for the implementation of the invention comprising, 
generally speaking, a control element I including a gestural 
sensor 1 controlled by an operator 10 and transforming the 

10 actions thereof into control parameters a> r , y, L t a 

modelization element II on which the control parameters act, 
including one non-linear part 2, associated with a linear part 3, 
and an element III creating the sound, including a means 4 for 
generating, from signals computed by the modelization element 

15 II, a signal which is transformed into sound synthesised by a 
digital/analogue converter 5. As known, the physical phenomena 
involved during the production of the sound of the clarinet, are 
expressed on the one hand, by a linear propagation equation of 
the waves in the pipe with loss, and on the other hand, by a 

20 non-linear equation linking the flow with the pressure and the 
displacement of the reed at the mouthpiece of the instrument. 

A simulation model of the sound therefore includes a 
linear part of the model corresponding to the resonator of the 
instrument which, in the case of the clarinet is composed of a 

25 cylindrical tube. For such geometry, assuming that the radius of 
the tube is large relative to the thickness of the boundary 
layers, the acoustic pressure inside the tube is governed by an 
equation in the form: 

3Q d 2 p(x,t) 1 d 2 p(x,t) a}p(M) =0 

dx 2 c 2 at 2 at | 



where 




R being the radius of the 



tube, i.e. 7mm in the case of the clarinet. The values of the 
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physical constants, in mKs units, are: c=340, lv=4.10-8, 

lt=5.6.10-8, Cv 

While looking for the solutions of the exp (i(cot-k(co)x)) 
type, co being the angular frequency of the wave, the following 
5 may be written: 



k(cD) 2 =^(l-iKac 2 GTX) (1) 



by replacing Vl + x with the approximated value 1 + f> 

10 when x is small, the conventional approximated expression of 
k(co), which will be used below, becomes: 



k(co)=: ^_^ acco ^ (2) 
c 2 

15 We know that, if we consider a pipe of infinite length, 

supposedly excited at x=0 and t=0 by a unit pulse 8(x)5(t), in all 
points x>0, the acoustic pressure propagated from such source 
is written in the form of a continuous sum of all the waves likely 
to propagate in the pipe: 



p(x,t) = Jexp(-ik(co)x)exp(icot)dco 



which appears as the inverse Fourier transform of the value 
exp(-ik(co)). 

25 In the following, the term "waveguide" will be reserved for 

the so-called Green formulation representing the propagation of 
a wave in a medium, and including the dissipation and the 
dispersion. 

In this Green formalism, the transfer function of a pipe of 
30 length L, representing the propagation, the dissipation and the 
dispersion is : 



F(co) = exp(-ik(co)L) = exp 




exp 




(3) 



10 



The dissipation, represented by the modulus of F(co), and 
the dispersion, represented by the phase of F(a>), are therefore 
proportional to Vco, while the propagation delay is provided by 

— . The length of the pipe will be therefore the control parameter 

c 

of the height and its radius the control parameter of the losses. 

We know, on the other hand, that the Fourier transform of 
the dimensionless pressures and flows at the input (Pe(oo), 
Ue(co)) and at the output (Ps(co), Us(co)) of the resonator are 
linked by the system of equations : 



Pe(co) = cos(k(co)L)P s (co)+ / sin(k((o)L)U s (co) 
U e (co) = /sin(k(co)L)P s (oo)+ cos(k((o)L)U s (co) 

15 

Conventionally, in order to model the internal acoustic 
pressure, the radiation may be neglected. The open end of the 
instrument is therefore perfectly reflecting, which involves that 
Ps(co)=0. This enables to express the relation between pressure 
20 and flow at the input of the resonator : 

Pe(co) = / tan(k(co)L)U e (co) = Z e (co)Ue(co) (4) 



where Ze(co) = i tan(k(co)L) is the normalized input 
25 impedance of the resonator. 

In the case of a conventional single-mode reed or lips 
model, the dimensionless displacement x(t) of the reed with 
respect to its point of equilibrium, and the acoustic pressure 
pe(t) at the origin thereof, are linked by the equation: 



30 



i d^ + *d£) +x(t) (t) 

mf dt 2 <°r « 



with the sign + when the pressure tends to close the reed 
or the lips and the sign - when the pressure tends to open it, 

and wherein co r =27i/ r corresponds to the resonance 
frequency /r, for instance 2500 Hz and q r is the quality factor of 
5 the reed, for instance 0.2. 

By writing, the equation (5) with the sign + in the Fourier 
domain, we obtain the transfer function of the reed: 

X(co) _ CD* 



Pe(co) G) r 2 — co 2 + icoq r co r 
10 whereof the angular frequency response is provided by: 

x(t) = 2( ° r expf- 1 co r q r tl sinf^ ^4 - qr 2 co r t | (7) 



V4 - qr 2 I 2 
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For exemplification purposes, Figure 2 provides two 
is diagrams respectively indicating, on the left, the transfer 
function and, on the right, the pulse response of the reed model, 
for a resonance frequency fx- 2 500Hz and a quality factor 

q r =0,2. 

As shown below, it should be noted that x(0)=0. 
20 Besides, in the case of a clarinet-type reed instrument or 

trumpet-type mouthpiece instrument, the acoustic pressure pe(t) 
and acoustic flow ue(t) (dimensionless) at the input of the 
resonator are linked in a non-linear manner by the equation: 

25 u e (t) = 1 (1 - sign(y - x(t) - 1))sign( Y - p.(t)K(1 - y + x(t)) A /|y - p e (t)| (8) 



In the case of a reed instrument, the parameter £ is 
characteristic of the mouthpiece and takes into account the 
position of the lips and the section ratio between the bill and the 
30 resonator. Such parameter £ is proportional to the square root 
of the opening of the reed in idle position and usually ranges 
between 0.2 and 0.6. 



20 
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The parameter y is the ratio between the pressure inside 
the mouth of a player and the plating static pressure of the 
reed. For a no loss pipe, it ranges from i for the initiation of 

vibrations to ^ for the position of a beating reed. 

5 The parameters £ and y are therefore two important 

playing parameters insofar as they represent, respectively, the 
way the player pinches the reed and the pressure of the breath 
into the instrument. 

While combining the displacement equation of the reed or 

io of the lips, the impedance relation and the non-linear 
characteristic, it appears that the acoustic pressure and flow, at 
the mouthpiece, are controlled by the following system of 
equations: 

1 d 2 x(t) q r dx(t) /4A ... /m 
15 — — + — — J 7 2 + x(t) = ±p e (t) (9) 

co r dr co r dt 



P e (co) = i tan 



a coo 

V C 2 



U e (co) (10) 



u e (t) = 1(1 - sign(y - x(t) - 1))sign(y - p e (t)K(1 - y + x(t))J|y - p e (t)| (11) 



The aim of the invention is therefore to find a formulation 
of the impedance relation in the time domain enabling to solve 
this three-equation system, by modelization the impedance 
relation in terms of elementary waveguides. 
25 To model the input impedance of the resonator in terms of 

elementary waveguides, the Fourier transform of the impedance 
Z e (co) is written in the form: 

Z fco) - itanrk(o>)L) - i sin ( k (°>> L > - exp(ik(co)L) - exp(-ik(co)L) 
Z.(a» - itan(k(a»L) " ' cos(k((o)L) - e xp(ik(co)L) + exp(-ik(co)L) 



30 




This expression may be written in the form: 

Pe(co) = 1 exp(-2ik(co)L) 

° } " Ue(co) 1 + exp(-2ik(co)L) 1 + exp(-2ik(o)L) 



5 Figure 3 represents a calculation model by combination of 

waveguides, directly derived from such last equation and whose 
transfer function is the input impedance of the resonator. It is 
composed of a sum of two elementary waveguides. The upper 
element corresponds to the first term of the equation (12) while 
10 the lower element corresponds to the second. The filter whose 
transfer function is ~ = - exp(-2ik(co)L) represents a two-way 
travel, with a sign change of the acoustic pressure at the open 
end. 

For exemplification purposes, Figure 4 provides two 

15 diagrams representing respectively, for a pipe of length L = 
0.5m and of radius R=7mm, at the top, the variation of the input 
impedance of the resonator in relation to the frequency and, at 
the bottom, the angular frequency response of the 
corresponding waveguide model, calculated by inverse Fourier 

20 transform of the impedance. 

The system of the three coupled physical equations (9), 
(10), (11) enables to introduce the non-linearity in the form of a 
loop connecting the output pe of the resonator to the input ue. 

Figure 5 provides an equivalent calculation diagram 

25 enabling, for the simulation of a reed instrument or mouthpiece 
instrument, to couple in a non-linear way the displacement of 
the reed or of the lips and the acoustic pressure with the 
acoustic flow at the input of the resonator, by computing, at 
each sampled time, the internal acoustic pressure at the 

30 mouthpiece. 

The model is entirely driven by the length L of the 
resonator and at least two parameters £ and y representative of 
the non-linear physical interaction between the source and the 
resonator, by means of a loop connecting the input to the output 




of the linear part and comprising a non-linear function playing 
the part of an excitation source for the resonator. As shown on 
Figure 4, the linear part takes up the diagram of Figure 2 and 
the non-linear function f is controlled by both parameters £ and 

5 y enabling to simulate the performance of a player, and has as 
input parameters, in the case of a clarinet, the pressure at the 
mouthpiece and the displacement x(t) of the reed with respect 
to its equilibrium point, computed in relation to the pressure at 
the mouthpiece, by a reed model (m) which forms the exciter. 

10 No input signal is necessary, insofar as y is directly 

proportional to the pressure in the mouth of the player. It is 
therefore the non-linearity itself, and its evolution imposed by 
the player, which plays the part of the excitation source, 
according to the physical model. 

15 For real-time synthesis of the sounds to be simulated, the 

model requires digital sampling and, to do so, a formulation is 
prepared, in the time domain, of the angular frequency response 
of the resonator, corresponding to the inverse Fourier transform 
of the impedance. Such formulation in the time domain enables 

20 to calculate the pressure pe(t) at the mouthpiece in relation to 
the flow ue(t) but, to do so, it is necessary to approximate the 
losses represented by the filter F(co) by means of an 
approximated digital filter. We shall therefore use an 
approximation of the transfer function of the filter, 

25 - F(co) 2 = -exp(-2ik(a>)L), where the coefficients are determined 

from physical variables such as the length of the resonator and 
its radius, in order to be able to make the necessary 
modifications in relation to the geometry of the resonator. In this 
view, we express analytically the coefficients of the digital filter 
30 as functions of physical parameters. 

In practice, it is particularly advantageous to use a single- 
pole filter while expressing the transfer function approximated in 
the form: 
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~ = b 0 exp(-2i CT D) 
. 1 - a, exp(-iro) v ' 

CO 

wherein ' e , fe being the sampling frequency, and 

D = f e — the pure delay corresponding to an away or return 
c js 

travel of the waves in the resonator. 
5 The parameters bO and a1 are expressed in relation to the 

|F(co) 2 | 2 = |F(uj)| 2 

physical parameters so that 1 1 1 1 for two provided 
values of co. 

The first value adopted col is that of the fundamental play 
frequency. This enables to ensure a down-slope time of the 
10 fundamental frequency of the angular frequency response of the 
waveguide model using the approximated filter, identical to that 
of the waveguide model using the exact filter. 

The second adopted value co2 is that of a harmonic 
selected in order to obtain global identical decrease in the 
15 angular frequency responses of the waveguides, respectively, 
exact and approximate waveguides. 

The selection of such second value co2 is therefore less 
restricted. It corresponds, for instance, to the second resonance 
peak in the case of the clarinet but, in certain cases, as shown 
20 below in the case of the trumpet, it may be preferable to select 
a higher rank harmonic. 

The system of equations to be solved is provided by: 
|FK) 2 | 2 (1 + a? -2a, cos^)) = b 2 



25 |F(co 2 ) 2 | 2 (1 + a 2 - 2a, cos(<u 2 )) = b 



wherein F(co) 2 = exp 



v v ^ J 

By neglecting the dispersion introduced by the non-linear 
part of the phase of F(co), the frequencies of the harmonics are: 
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where 2k-1 is the rank of the harmonic. 



Lfe 

By noting : 

i |2 i |2 

c, = cosCct^, c 2 = cosCtBj),^ = pfco,) 2 ! ,F 2 = |F(o) 2 ) 2 | , A, = F,c tl A 2 
the coefficients a1 and bO are provided by : 



a, = 



_ A1- A2- A /(A 1 - A 2 ) 2 -(F, -F 2 ) 2 
F, -F 2 



(14) 



b 0 = 



^2^(0, - c 2 )(A 1 - A 2 - ^/(A, - A 2 ) 2 - (F, - FJ) 

Ft -F 2 



(15) 



From the expression of the input impedance of the 
resonator: 



Ze(co) = 



1 



exp(-2ik(co)L) 

1 + exp(-2ik(co)L) 1 + exp(-2ik(co)L) 



and by noting z=exp(ixn), we obtain directly: 



P.(z) 
U e (z) 



p-2D 



1 - a^z 



-1 



1 + 



-2D 



1 + 



-2D 



{1 -a.z" 



1 - a,z- 1 - b 0 z- 2D 
1 - a.z" 1 + b 0 z' 2D 



wherefrom the differential equation is derived: 

Pe(n) = u e (n) - a lUe (n - 1) - b 0 u e (n - 2D) + a lPe (n - 1) - b oPe (n - 2D) (16) 
Figure 6 provides two diagrams respectively showing, at 
the top the variation of the input impedance of the resonator 
approximated in relation to the frequency and, at the bottom, 
the angular frequency response of the corresponding waveguide 



model, calculated from the differential equation for a cylindrical 
pipe having a length L=0.5m and a radius R=7mm. 

It can be noted that the model thus established provides a 
angular frequency response very close to that of the resonator, 
5 represented in Figure 2. 

As well as for the filter representing the losses, the 
relation between the acoustic pressure and the displacement of 
the exciter (reed or lips) must be discretized in the time domain. 
Or, the angular frequency response of the reed is a sinus wave 
10 damped exponentially which meets the condition x(0)=0, as 
already mentioned above. II is therefore possible to build a 
digital filter for which the displacement of the reed at the time 

t n = -p- is function of the acoustic pressure at the time 
fe 

n 1 

t n-1 = and not t n . This enables to verify the property x(0)=0 

fe 

15 of the continuous system when the reed is subjected to a Dirac 

excitation. To comply with such condition, instead of using the 

bilinear transformation to approximate the terms ico and -co 2 , we 

fe 

use, according to the invention, the expressions ico « — (z-z~ 1 ) 

and -co 2 ~f e 2 (z - 2 + z~ 1 ), which correspond to an exact centred 

20 digital second order differentiation. 

With these approximations, the digital transfer function of 
the reed is provided by: 



X(z) 
Pe(z) 





co r 2 +f e 2 (z 


-2 + z" 1 ) + — 
2 


(z - z" 1 )q r co r 






z" 1 




f e 2 + f e q r 


z" 1 


f 2f 2 "I 


z - 2 ff e q r U 2 ) 


co 2 2co r 




I W r J 


[2<o r co 2 J 



wherefrom the differential equation is derived: 
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x(n) =0xp e (n)+biaPe(n-1)+ai a x(n-1)+a 2a x(n-2) (17) 
wherein the coefficients b1a, a1a and a2a are defined by 



2f e 2 „ fq f : 

f e 2 fe<3r t 1 co r 2 2co r CD? 

a 0a = + ^ b 1a = — a 1a =— ^ a 2a = — <- l 

co r 2co r a 0a a 0a a 0a 



Figure 7 provides, for such a reed model approximated, 
two diagrams representing, on the left the transfer function and 
on the right the angular frequency response, with a sampling 
10 frequency fe=44 100 Hertz, the values of the parameters being 
fr=2500Hz and qr=0,2. 

It can be seen that the diagrams obtained are very close 
to those of Figure 2 

This being established, we shall now expose an explicit 
15 resolution method, according to the invention, for coupling the 
differential equations with the non-linear characteristic. 

The sampled formulations of the angular frequency 
responses of the displacement of the reed and of the impedance 
enable, indeed, to write the sampled equivalent of the system of 
20 equations (9, 10, 11) stated above, in the form: 

x(n)=bi a Pe(n-1 )+ai a x(n-1 )+a 2a x(n-2) (18) 

pe(n)=u e (n)-a 1 Ue(n-1)-boUe(n-2D)+aiPe(n-1)-boPe(n-2D) (19) 

25 

u e (n) = ±(1 - sign(y - x(n) - 1))sign(y - p e (n)K(1 - y + x(n))^/|y -p e (n)| (20) 

This system of equations is implicit, because the 
calculation of p e (n) by the impedance equation requires to know 
u e (n) and such value is itself obtained from the non-linear 
30 equation and requires to know p e (n). Still, as already mentioned 
above, the calculation of x(n) does not requires to know p e (n) 
but to know p e (n-1) which is indeed known at the time n. 
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This enables, according to the invention, to solve simply 
and exactly the coupled system. To do so, the terms of the 
equations 19 and 20 above, which do not depend on the time 
sample n, may, indeed, be gathered in the expressions: 

V=-aiU e (n-1)-boUe(n-2D)+aiPe(n-1)-b 0 Pe(n-2D) 



To generalise the method, it is interesting to associate u e 
(n) with a coefficient bc 0 = 1 in the case of a cylindrical pipe. 

Both equations 19 and 20 above may then be written in 
the form: 

Pe(n)=boC 0 u e (n)+V 

u e (n) = Wsign(y - P e (n))Jy - P e (n)| 

1 

Since the term - (1 - sign(y - x(n) - 1)) cancels W when (1- 

y+x(n)) is negative, W always remains positive. If we 
successively consider both cases: y-pe(n)>0 and y-pe(n)<0 
corresponding respectively to the cases ue(n)>0 and ue(n)<0, 
ue(n) may be expressed exactly and without involving the 
unknown pe(n), in the form: 



Consequently, the calculation of the acoustic pressure and 
the flow at the mouthpiece at a sampled time n, may be 
conducted by using sequentially the following equations: 



W = - (1 - sign(y - x(n) - 1)) x «1 - y + x(n)) 



which enables to write: 



p e (n) =u e (n)+V 




x(n)=biap e (n-1 )+a 1a x(n-1 )+a 2a x(n-2) 



(21) 



V=-aiU e (n-1)-boUe(n-2D)+aiPe(n-1)-boPe(n-2D) 



(22) 



W = -L (1 _ sign(y - x(n) - 1)K(1 - y + x(n)) 



(23) 




(24) 



p e (n)=bc 0 u e (n)+V 



(25) 



Thus, the invention enables to solve in the time domain, 
the system of equations governing the physical modelization of 
the instrument, from a sampled formulation equivalent to the 
angular frequency response of the displacement of the reed, of 
the impedance relation and of the non-linear characteristic, 
which is translated into the system of equations 18, 19, 20, 
wherein : 

- the equation (18) is a digital transcription of the model 
(m) of Figure 5, 

- the equation (19) is a digital transcription of the 
impedance model of Figure 3, 

- the equation (20) is a digital transcription of the non- 
linear characteristic linking the displacement of the 
reed and the acoustic pressure with the acoustic flow. 

By grouping the terms which do not depend on the time 
sample n, the method according to the invention enables, 
indeed, to determine the flow and the pressure at the input of 
the resonator by a sequential calculation of the equations 21 to 
25, and to solve, in the time domain, the system of equations 9, 
10, 11 governing the physical modelization of a clarinet-type 
reed instrument, in order to synthesise the sounds produced by 
such an instrument. 

For exemplification purposes, Figure 7 shows the variation 
of the internal acoustic pressure at the mouthpiece, calculated 




by such a non-linear model involving waveguides, for a pipe 
having a length L=0.5m and a radius R=7mm, the values of the 
parameters being 7=0.4, <;=0.4, fr=2205 Hz, qr=0.3. 

Three phases can be observed: the attack transient 
5 corresponding to sudden increase of y and the steady state 
while y and £ diminish gradually, in a linear manner, up to the 
oscillation threshold, and the extinction transient. 

In practice, the digital implementation of such a non-linear 
waveguide model, may be conducted with the use of elements 

10 available on the market for the gestural sensor. For instance, 
digital implementation is possible in language C in the form of 
an external « clarinet » object for the environment known under 
the trade name Max-MSP, driven from MIDI controls supplied by 
a controller Yamaha WX5®. This controller measures the 

15 pressure of the lips on the reed, which controls the parameter <;, 
and the pressure of the breath, which controls the parameter y. 
Such information received in MIDI format (therefore between 0 
and 127) are re-standardised to correspond to the scale of the 
physical parameters. The waveguide is tuned from the 

20 information MIDI pitch controlled from the finger position which 
determines the length L of the pipe. 

Still, as in a real instrument, the height changes in relation 
to the physical parameters such as y, cor and qr. Still, in 
reality, a musical instrument is not tuned perfectly for all the 

25 finger positions. The use of an all-pass filter for the 
implementation of the fractional part of the delay D is therefore 
not compulsory. 

In practice, it appears that the playing sensations of such 
a virtual instrument are quite comparable to those of a real 

30 instrument. 

Still, the method which has just been described for the 
simulation of the sound of a clarinet-type reed instrument, may 
still be perfected. 

We know, indeed, that the acoustic pressure at the 

35 mouthpiece is not the variable representative of the sound 




perceived. It is therefore interesting to calculate the external 
pressure which, for a cylindrical pipe, may be expressed as the 

time derivate of the outgoing flow : p ex t(t)= dUs ^ , By neglecting 

dt 

still the radiation, which involves p s (t)=0, it becomes: 

5 

Pe(co) = /sin(k(co)L))U s (co) 
Ue(co) = cos(k(co)L)U s (co) 
wherefrom we may derive: 

10 Us(co) =exp(-ik(co)L)(Pe(co)+U e (co)) 

From the perception viewpoint, the term exp(-ik(co)L) is 
negligible. The expression above may therefore be simplified 
and becomes: 

15 

Pext(t)=^(p e (t) + u e (t)) (26) 
dt 

Thus, from a digital viewpoint, the calculation, at each 
sampled time n, of the external pressure p ex t( n ). is reduced to a 

20 simple differential between the sum of the internal pressure and 
of the flow, between the time n and the time n-1. 

As shown in Figure 1, which represents schematically a 
whole digital instrument for the implementation of the invention 
in the case of a wind instrument, the signals p e (t) and u e (t) 

25 enabling the calculation of the external pressure p ex t(t) are 
prepared by the modelization element II from control parameters 
CDr, C> y, L. 

For a cylindrical resonator, such modelization element II is 
of the type represented in Figure 5 and enables coupling the 
30 three equations (9), (10), (11). 

The linear part 3 includes a computing block 31 of the type 
represented in Figure 3, the transfer function Ze(co) of which is 
the input impedance of the resonator. 



The model is driven by the length L of the resonator and 
the non-linear part 2 implements a non-linear function 21 
controlled by both parameters £ and y and having as input 
parameters the pressure p e (t) calculated by the linear part 3 and 

5 the displacement x(t) of the exciter 22 calculated, in the case of 
the clarinet by a reed model (m) in relation to the same 
pressure p e (t) at the mouthpiece. 

From such pressure p e (t) and the flow u e (t) at the 
mouthpiece, computed respectively by the linear part 3 and the 

10 non-linear part 2, the bloc 4 computes the sound signal p ex t(t) 
emitted by the digital instrument thanks to the converter 5. 

The method according to the invention has been described 
in detail for the simulation of a simple instrument, with 
cylindrical resonator, of the clarinet type. 

is Still, as will now be seen, the general diagram of Figure 1 

may be applied to the simulation of more complex phenomena. 

In particular, physical measurements have shown that the 
vibrations of a reed are more complex than a simple sine wave 
damped. It has therefore appeared that the simple model 

20 described above could be perfected in order to improve the 
quality of the sound produced, as perceived. To do so, we 
consider a very simplified reed model in the form of a free 
embedded string. Thus, the impedance model of a cylindrical 
pipe described above, which generates mainly odd harmonics, 

25 acts as a basis for the realisation of a multiple mode reed model 
complying with the same condition x(0)=0, which enables to 
keep the digital resolution diagram of the equations 21 to 25 
specified above. 

Since the value of the impedance is real for all the 

30 resonance peaks, its angular frequency response is a sum of 
cosine functions. On the other hand, it has been seen that the 
angular frequency response of a single mode reed model is a 
sine wave function. We shall therefore use, in the definition of 
the model, the transfer function between a damped sine wave 

35 and a damped cosine wave. While keeping the same notations, 
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the Fourier transform of x(t)=-;==exp 



V4-qr 2 



-^a»rqr* ) sin [^ V 4 " 9r 2<B rt j 



X(o) = — 
is provided by : ^ " 03 + ,coq ^ 



Similarly, the Fourier transform of 
/, 2t ° r , ex P| ~ ^ ^rPrt ] cos 



1 V^^m] 



is provided by : 



Y(a>) = 



co r (co r q r +2ico) 



^4 - q r 2 (co r 2 - <B 2 + ia>q r co r ) 



The transfer function between y(t) and x(t) is then: 
X(co) = J 4 - q r 2 co r 
Y(co) co r q r + 2ico 

io We may thus write the transmittance model in the form: 



X(co) 
P.(<°) 



c 


1 - B exp 


r Y 

. 7100 


) 




+ ico 

J 


1 + B exp 

V 




) 



To determine the three unknown variables coa, C, B, three 
is conditions are set. The first condition consists in keeping the 
frequency of the first peak, which is the maximum peak. To do 
so, we select a>a = cor assuming that the frequency offset 
resulting from the quality factor qr can be negligible. 

The second condition, satisfied by the single mode reed 
20 model, consists in imposing a unit value of the transmittance 
module for co=0 in order to maintain X(co) * Pe(co) at low 
frequencies. 

The third condition is an imposed value of — for the 
transmittance module at the frequency cor, in order t9 r keep the 
25 height of the peak of the single mode reed model. 
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Thanks to these three conditions, the transfer function 
thus realised reproduces the main characteristics of the single 
mode reed model. 

Both first and second conditions lead to the system of 
5 equations: 

2C(1 + B) 1 



CO 



rA /q r 2 +4(1-B) q r 
2C(1 - B) 



= 1 



<B,q,(1 + B) 

2 1 

assuming: A, = A 2 =— co r q r and A3=A 2 A 1 q r , the 

«Wqr 2 + 4 2 

coefficients B and C solving the system are provided by: 

B= 1 + A 3 -2VA7 
1-A, 

c= Aa_~ VAT (28) 



(VAT - l)qA 



Similarly to the case of the single mode reed, we shall 
prepare a digital model wherein the displacement of the reed at 

20 the sampled time t n =^- is a function of the acoustic pressure, not 

n 1 

at time tn but at time t n -i= — — . This is possible since the angular 

frequency response of the model is the sum of damped sine 
wave functions. 

To fulfil such condition, we use an approximation of ico in 
25 the form ico * f *( z " 1 > = f e(exp(\w) - 1) 

Moreover, to add to the model an additional control 
parameter, we replace the coefficient B with a filter such as 



10 



15 
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B« - . Thus, we may adjust the damping of the 

1 - a a exp(-im) 

harmonics in relation to the damping of the fundamental. To 
keep the characteristic X(0)=1, the parameters b a and a a are 

linked by the equation B= a 



1-a a 



The term expf-i — ] is replaced with its sampled 
equivalent: z" Da =exp(-ituD a ) with the delay D a defined by: 



D a =E 



— j wherein E indicates the integer part. 



1 

Assuming: p = -oo r q r , the digital transfer function is 
written: 

X(z) = C(z Da+1 - a a z Da - b a z) 

P e (z) (p + fe(z - 1))(z Da+1 - a a z Da + b a z) 

= C(z- 1 -a a z- 2 -b a z- Da ' 1 ) 

f. - (f.(1 + a.) - p)z~ 1 -a a (p - f e )z~ 2 + f e b a z' Da - b a (f e - P)z- Da ~ 1 

which leads to the differential equation : 

x(n)=b a ipe(n-1 )+ba2Pe(n-2)+baDip e (n-D a -1 ) +a a ix(n-1 ) +aa2x(n-2)+aaDx(n-D a )+aaDix(n-Da-1 ) (29) 

wherein the coefficients a a i, a a2 , a aD 2, a aD i are defined by: 

a „= ifliiLJ , ^^ilA . a>D= . b „ to =ilki 

'e 'e 'e 

and the coefficients b a1l b a2 , b aD i by: 

ba1= Y~ ' b a2= — — b a D1 = — . 

f e fe fe 
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We notice that equation 29, thus established, enables to 
determine the dimensionless displacement x(n) of the reed at 
the sampled time n, from previous times. 



It appears therefore that, in the case of a multimode reed 
model, the digital calculation diagram of the sound may be the 
same as that of the single mode reed while replacing equation 
21 with equation 29 which is another digital transcription of the 
5 model (m) of Figure 5. 

For exemplification purposes, Figure 9, which is 
analogous to Figures 2 and 7, provides, for an approximated 
multimode reed model computed according to the invention, two 
diagrams representing respectively, in solid lines, on the left, 

10 the transfer function and, on the right, the pulse response with a 
sampling frequency f e =44100 Hertz, the parameters having 
values f r =1837,5 Hz, q r =0,2, a a =0. 

On the same diagrams, we have superimposed, as doted 
lines, the transfer function and the angular frequency response 

is of the single reed model as represented in Figure 7. 

According to another development of the invention, we 
may also improve the sound model of clarinet in order to make it 
more natural by integrating a certain noise thereto, the system 
being thus more realistic. Since the noise is created by a 

20 turbulence at the reed before the beginning of the pipe, the 
noise is added to x(t). It appears on the other hand that, in 
practice, the level of noise depends on the pressure of the 
breath whereas its « colour » depends on the pressure of the 
lips on the reed. Indeed, from a physical viewpoint, the harder 

25 the reed is pressed, the smaller the opening between the reed 
and the pipe and the greater the turbulence. Therefore we shall 
use a simple noise model whose level is driven by y and the 
brilliance driven by £. 

For the digital implementation of such a noise model, we 

30 shall use a low-pass filter of a white noise. The transfer function 

of this filter is provided by T b (z) = bb ^ 1 Z a »\ the coefficient bb 

1 - a b z-' 

being driven by y, and the coefficient ab driven by The 
variation laws of bb and aa may be determined so that the 
sound simulated by the model is as realistic as possible. 




Both diagrams of Figure 10 show for exemplification 
purpose, the variation of the module of the spectrum of the 
external acoustic pressure corresponding to the sound produced 
by the model, respectively on the top diagram for a single mode 
5 reed and on the bottom diagram for a multiple mode reed with 
additional noise, the simulation parameters being as follows: 

f r =2205 Hz, q r =0.25, a a =0, y=0.44, ^=0.4, L=0,48 m, R=7.10" 3 m 

10 As already mentioned, the method according to the 

invention, as it has just been described in details, refers to the 
simulation of sounds produced by a reed and cylindrical 
resonator musical instrument of the clarinet type. But the 
invention is not limited to such an application and may, 

15 conversely, be subject to numerous developments. 

Indeed, from the non-linear physical model, involving 
waveguides and schematised in Figure 5, as well as of its digital 
transcription according to the sequence of equations 21 to 25, it 
will be possible, to simulate the operation of a resonator having 

20 any geometry, by modifying the impedance model for cylindrical 
resonator and the associated differential equations, to replace 
the latter with impedances and more complex differential 
equations. 

While keeping the properties of the model which has just 
25 been described, it is possible indeed, from the general diagram 
of Figure 1, to build other complex impedance models, by 
combining certain impedance elements in parallel or in series 
and by offering digital approximations enabling explicit usage of 
the physical variables and more flexible control of the digital 
30 instrument. 

For exemplification purposes, we shall now describe 
certain developments of the basic model for cylindrical 
resonator, with reference to Figures 11 to 14 which represent 
equivalent calculation diagrams involving waveguides and 
35 corresponding to resonators having diverse geometries. 
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Generally, on these diagrams which correspond to the 
computing block 31 in Figure 1, the operator C(co) represents 
the input impedance and C' 1 (oo) the input admittance of a 
cylindrical resonator, the digital model corresponding to C" 1 (co) 
5 being obtained by changing the sign of the coefficient bo only. 

A first improvement on the basic model which has just 
been described with reference to Figures 3 and 5, will enable, 
by the use of wave guides similarly, to realise a physical model 
for cylindrical resonator with terminal impedance. Such an 
10 element will enable, for instance, to link together parts of 
cylindrical resonators having different lengths and sections, in 
order to simulate the input impedance of a conduit of variable 
section, or still to take into account the radiation impedance. 

To do so, we consider the formalism of the transmission 
15 line linking the acoustic pressure and flow, respectively to the 
input of the resonator (Pe(co), U e (co)) and to its open end (P s (co), 
Us(co)). 

Noting Z c = the impedance characteristic, we have: 

20 P e (co) = cos(k(co)L)P s (co) + iZ c sin(k(co)L)U s (co) 

Ue(co) = -L sin(k(co)L)P s (co) + cos(k(co)L)U s (co) 

By noting Z s = the output impedance and 

U s (co) 

25 R(co) = — — , we can write the input impedance in two 
' Z c + Z s (co) 

different ways: 

Z s (co) 



+ itan(k(co)L) 

P -^> - Z < (30) 



Z c U e (a>) 1 + j Z>) tan(k(to)L) 




= 1 - R(co) exp(-2ik(co)L) 
1 + R(co) exp(-2ik(co)L) 

The equation 31 shows therefore that the impedance of a 

cylindrical resonator with terminal impedance may be obtained 

from the impedance of a cylindrical resonator without terminal 

5 impedance, while replacing: 

exp(-2ik(co)L) with R(co)exp(-2ik(co)L). 

Figure 11 provides an equivalent calculation diagram 

involving waveguides, for the implementation of the equation 

30, enabling to calculate the impedance of a cylindrical 

10 resonator with terminal impedance. 

Such a model enables to generate in cascade the input 

impedance of a conduit having any geometry and liable to be 

defined by a succession of elementary cylindrical conduits. 

Consequently, the invention may be applied to the 

15 simulation of the vocal conduit. 

But the invention may be subject to many other 

developments and the diagram of Figure 1 enables in particular, 

from the basic physical model for cylindrical resonator 

schematised on Figure 5, to build specific models for the 

20 simulation of diverse musical instruments. 

It is thus that in a first development of the basic model, we 

shall build a model for conical resonator, usable, for instance, 

for the simulation of a saxophone. 

If we designate by L the length of the pipe, R its input 

25 radius, 9 its opening, the distance x e between the apex of the 

cone and the input is: 
R 

The input impedance, relative to the characteristic 
30 impedance 

7iR 2 j S then provided by the expression: 
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Pe(C0) _ 1 _ (32) 



Z c U e (co) 1 , 1 

i tan(k(co)L) 



I co — 
c 



which may be written in the form: 
x e 

D / \ ICO — 

P e(<°) _ C 



Z c U e (co) ^ 
1 + 



i tan(k(co)L) 

Figure 12 provides an equivalent calculation diagram 
involving waveguides, wherein the element noted D is the 
differential operator D=ico and the element noted C-1(co) 
corresponds to the diagram of Figure 2 by replacing 
-exp (-2ik(co)L) by +exp (-2ik(co)L). 

The corresponding differential equations may be 
established by following a process similar to that which has 
been described above. 

Thus, using the bilinear transform to approximate ico, the 
digital transfer function of a conical resonator is provided by: 
P.(z) 1 



Z c U e (z) z + 1 1-a^" 1 +b 0 z 



-2D 



2fe~(z-1> 1 " a,z" 1 - b 0 z- 
c 

1 1 

By noting : Gp = 1+ — — and Gm = 1- — — , the 

2f e ^ 2f e ^ 

c e c 

transfer function is reduced to: 

P.(z) 1 - (a 1 + 1)z- 1 + 3l z" 2 - b 0 z" 2D + b 0 z- 2D " 1 



Z c U e (z) G p - ( 3l G p + G m )z- 1 + a.G^" 2 + b 0 G m z" 2D - b 0 G p z 
wherefrom the differential equation is derived: 



-2D-1 



p e (n) = bCoUe(n) +bCiUe(n-1)+bC2Ue(n-2)+bc D Ue(n-2D)+bCDiUe(n-2D-1) 
+ ac 1 p e (n-1)+ac2Pe(n-2)+ac D Pe(n-2D)+acDiPe(n-2D-1) (33) 
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wherein the coefficients bcO, bc1, bc2, bcD and bcD1 are 
defined by: 

1 a i + 1 3i b 0 b 0 

5 bc 0 = — , bci=- , bc 2 =-^, bc D =- bc D1 =-^ 

G p G p G p G p G p 

and the coefficients ad, ac2, acD and acD1 are defined 

by: 

ac 1 — , ac 2 - — — , ac D - — — , ac D1 - d 0 

G p G p G p 

10 Similarly, the invention may be applied to the case of 

short resonators which appear, for instance, in the mouthpiece 
of a brass instrument or in the bill of a reed instrument, or of a 
register hole or lateral hole. 

To do so, it is admitted that the radius of the short 

15 resonator is sufficient enough to keep the loss model used until 
now. 

This approximation of a short resonator will consist in 
approximating the impedance ZI(co)= i tan (k(oo)l) for small 
values of k(co)l. 
20 We thus obtain the expression 

Z,(co)= i tan (k(co)l) = l(k(o)l) = G(co) + icoH(co) (34) 

and H (©)=-(l- G(co)). 



1 - expf- ac^l 
wherein G(co) = 



1 + exp 



-ac, 

v 

25 In another development of the basic model for cylindrical 

resonator, the invention also enables to simulate a more 
complex resonator, by assembling elementary impedances 
representing, in the one hand, the conduit and, on the other 
hand, the bill of a reed instrument or the mouthpiece of a brass 

30 instrument. 



In such a case, we shall model the mouthpiece or the bill 
by a Helmholtz resonator comprising a hemispheric cavity 
coupled with a short cylindrical pipe and a main resonator with 
conical pipe. 

The input impedance of the resonator assembly may be 
expressed by: 



Ze(co) = 



1 



. V 

ico + - 

pc 2 



ico- 



xe 



Z 1 (k 1 (co)L 1 ) + Z 2 



ICO - 



xe 



1 + 



itan(k2(co)l_2) 

wherein V=— tiR* is the volume of the hemispheric cavity, 
6 

L1 is the length of the short pipe, L2 is the length of the conical 
pipe, Z1 and Z2 are the characteristic impedances of both pipes 
which depend on their radii, k1(co) and k2(co) take into account 
the losses and of the radius R1 and R2 of each pipe. 
This equation may be written in the form: 



ico 



iZ l (k l (a))L l ) + Z 2 



ico - 



1+ 



itan(k 2 (co)L 2 ) 



\ + ia> 



10) ■ 



iZ ] (k l (co)L 1 ) + Z 2 



1 + - 



c 



/tan(A: 2 (6;)L 2 ) j 



We may thus establish the equivalent calculation diagram 
represented in Figure 13, wherein the operator noted Ci(co) 
corresponds to the diagram in Figure 5 and the operator noted 
S(co) corresponds to the input impedance of the conical pipe and 
to the diagram of Figure 12. 



20 



25 
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Resolution methods similar to those which have been 
described previously, enable to express the equivalent digital 
model and the corresponding differential equations. The 
impedance of the shoct pipe may be modelled by using the 
5 approximation expressed by the previous equation (34). The 
impedance of the conical pipe is represented by the model 
corresponding to the differential equation (33). The admittance 

of the cavity ico-^- is approximated by the bilinear 

pc 2 

z - 1 

transformation d , where d=2f e . 

z + 1 

io Considering the association of the hemispheric cavity and 

of the short pipe as a Helmholtz resonator having as a 

resonance frequency <o h = ■ we ma y use such frequency co h 



to approximate G(co) by G(co h ) and H(a>) by H(« h )- Both 
frequencies used for the calculation of the coefficients ai and b 0 

c(12ttL + 9n 2 x e + 16L) . 4 

15 are co, = — ~ which correspond to the first 

1 4L(4L + 3nx e + 4x e ) 

impedance peak of the conical pipe, and co 2 = co h - moreover, we 

shall use z n = — = Z 2 to normalize the input impedance. 

S 2 

The digital impedance of a brass-type resonator is thus 
provided by the expression: 



*■» ■ Vd(z-1) + ~ 1 (35) 



pc 2 (z + 1) ^ Ci(z) + ^s 2(z) 

^1 ^2 



which may be simplified in the form 

,k=4 be 

Ze(z) = ^ ° 



Ztl bc k z- k + Z^lg bc Dk z- 2D - k 
ac 0 — Zk^ 4 ac k z k - Zk = o ac Dk z 



wherefrom the differential equation may be derived : 




k=4 k=3 

p e (n)=£bc k u e (n - k) + £bc Dk u e (n - k - 2D) 

k=0 k=0 

+ £ac k p.(n - k) + X ac DkPe (n - k - 2D) (36) 

k=1 k=0 

5 

The coefficients are derived from direct calculation of the 
equation (35). 

The invention may still be applied to the modelization of a 
10 cylindrical resonator with register holes. To that effect we use 
elements involving waveguides and corresponding respectively 
to a physical model of a cylindrical pipe with terminal 
impedance representing a pipe of length L1 between the 
mouthpiece and the register hole, a model of short pipe which 
15 represents the register hole of length h t and the basic model for 
cylindrical pipe representing a pipe of length L2 between the 
register hole and the open end. 

The terminal impedance of the first part of the pipe may be 
written: 

1 

20 Z S (C0)= - 

Z t C t (<o) + Z 2 C 2 (co) 



By writing such expression in the form : 



<-2_ ^ 2 \™) 
Z t C» 

25 we may establish a general calculation diagram by 

waveguides represented on Figure 14, which, in the case which 
has just been described, is used to calculate the terminal 
impedance by a parallel combination of the impedances Z 2 C2(co) 
and ZtCt(co). 
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In the limit case of a closed register hole, with two parts of 
the pipe having the same radius, we assume Z t Ct(co)=oo 
and Z 2 =Z C 

leading to Z s (a)) = iZ c tan(k(co)L 2 ) and Z c (co) = iZ c 
5 tan(k((o)(Li + L 2 )). 

The total input impedance of the pipe may then be 
expressed by : 

P e (p) = (C 2 (co) + C>))Z 1 C i (Q) + Z.C.MC^co) 
Z c U e (co) (1 + C>)C 2 (a>))Z t C t (co) + Z c C 2 (co) 

Using ico=f e (1 -Z~ 1 ) in the impedance of the short pipe, the 
differential equation can be derived directly as previously. 

Thus, the model for simulating the cylindrical resonator of 
is the clarinet-type, obtained by direct transposition of the 
simplified equations of the physical behaviour of the instrument, 
may be adapted to the simulation of instruments with non- 
cylindrical resonator, such as saxophone, trumpet or other wind 
instruments. 

20 But the invention is not limited to such an embodiment and 

to adaptations which have just been described since, without 
departing from the scope of the claims, it may be applied to the 
simulation of other types of instruments, for instance when a 
string is rubbed as with the violin or struck as with the piano. 

25 Indeed, in the case of a string, it is known that using the 

formalism of mechanical transmission line, analogous to that of 
an acoustic transmission line, we establish analogous relations 
between the variables of the impedance relation. Said relations 
are the strength exerted on the mechanical element, and the 

30 speed of such an element resulting from such strength. Insofar 
as the quantity linked to the sound emitted is linked with the 
speed (i.e. the effect) in the mechanical case, by opposition to 
the acoustic case wherein the quantity linked with the sound 



5 



10 
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emitted is the pressure (i.e. the cause), we prefer to describe 
the resonator in terms of admittance rather than impedance. 

The equations of the mechanical transmission line 
between a point (b) and a point (e) are then: 



,F;(^0 = cos(/i:(u;)./:)F 6 (w) + iZ c sin (k(uj) ,L)V b {^) 



(38) 



wherein F and V represent respectively the strengths and 

speeds at each point. 

The wave number k(oo) is expressed conventionally from 
15 the differential equation of the movement of a string under 

deflection and comprises, as in the acoustic case, propagation 

(delay), dissipation, dispersion parts (see for instance : C. 

Valette, C. Cuesta "Mecanique de la corde vibrante", Hermes, 

treatise on new technologies, serie Mecanique. 1993). 
20 If we assume that the end (b) of the string is fixed and the 

end (e) mobile, we have: 



V b (u) = 0 V;M = ^4 = i-*tan(*(u;)L) 



This relation forms the input admittance of a part of 
embedded-free string at the point where it is free, and is 
identical, within one multiplying constant, to the acoustic 

30 impedance of a cylindrical resonator. It may therefore be 
represented by a diagram analogous to that of Figure 3. 

At the point where the interaction with the exciter is 
realised, we write the equations of continuity between both 
parts of string 1 and 2, by considering that the total strength 

35 exerted is the sum of the strengths exerted on each part, while 
the speeds of each part are equal: 



F(cj) = F 1 (lj) + F 2 (v) 
40 v(cj) = V x (cj) = V 2 (w) 




Assuming that the complete string is fixed at both ends, 
this enables to express the input admittance of the string from a 
series combination of each part of the string, i.e., in terms of 
admittance: 



10 

Such a relation is identical, within one constant, to the 
parallel combination of two cylindrical acoustic resonators and 
may therefore be represented by the diagram of Figure 14 
wherein the elements noted Ct and C2 represent the admittance 

15 of each part of the string. 

For exemplification purposes, Figure 15 represents, in 
relation to the frequency, at the top, the exact admittance of a 
string at the eighth of its length, calculated with an expression 
of k(a>) derived from a conventional model, and, at the bottom, 

20 the admittance approximated using an approximation of the 
losses with a digital first order filter whose coefficients are 
computed with the same method as in the acoustic case. 

It should be noted that, in the case of the violin, the 
contact point between the bow and the string is very close to 

25 one of the ends of the string. Consequently, it is possible to use 
the « short pipe » approximation for the admittance of either of 
both parts. Moreover, insofar as the losses expressed by k(co) 
are very small in a string, it is also possible to neglect it for the 
short part. Thus, the admittance of a violin string at the contact 

30 point with the bow may be expressed in the same manner as the 
impedance of an acoustic conical resonator and can therefore 
be represented by the diagram of Figure 12. 

As in the case of acoustic self-oscillating instruments, the 
admittance described in this basic model comprising a string 

35 with two fixed ends, may be refined in order to take into account 
additional physical phenomena. The method consists again in 
associating the admittances of different elements. Thus, it is 




possible with such an approach to build an input admittance of 
a resonator composed of two strings coupled by a harmonic 
table, i.e. whose the ends are not fixed any longer, but mobile. 
Indeed, in most notes of a piano, there are two or three strings 
5 tuned to very close frequencies, which are struck 
simultaneously by the same hammer. In such a case, the total 
admittance is expressed by an association of two admittances 
of identical strings, each of these admittances being composed 
of two parts of strings, one of these parts is expressed 

10 identically to the input impedance of a cylindrical pipe with 
terminal impedance. In the mechanical case, the terminal 
admittance corresponding to that of the sound board, may be 
expressed by combinations of localised elements similar to 
those employed to describe the mouthpiece or the bill (i.e. 

15 masses, springs, dampers), enabling to take into account one or 
several vibration modes of the sound board. 

According to the invention, the formulation of the 
resonator en terms of mechanical admittance may be 
processed, for instance in an instrument where the string is 

20 struck as with the piano. In such a case, and as in the acoustic 
case, the speed of a string struck by a hammer such as that of 
a piano, may be expressed from the system with three coupled 
equations: 



fit) = K(y h (t) - y s (t)ni - ct-(y h {t) - ! Js (t))) 




(40) 



V s (cj) = Y(u>)F(u>) 



35 i.e.: 

- a non-linear characteristic expressing the strength in relation 
to the relative displacements and speeds of the hammer and of 
the string, 




- an equation of the dynamics of the hammer linking its 
acceleration to the strength exerted by the string thereon by 
reaction, which is analogous to the expression of the 
displacement of the reed in relation to pressure, 

5 - an admittance equation expressing the speed of the string in 
relation to the strength which is imposed thereupon, which is 
equivalent to the acoustic impedance relation. 

The non-linear impact characteristic used here is known 
as the Hunt-Crossley characteristic. The exponent (p) ranges 

10 conventionally between 2 and 3, and is not an integer. yh(n) 
designates the displacement of the hammer, ys(n) that of the 
string. It should be noted that it is here a new way of writing the 
problem. Indeed, conventionally, the impedance relation used 
here is replaced with the differential equation of the movement 

is of the string. 

Such a simulation method of a string instrument may be 
implemented by a digital instrument whereof the model is 
schematised on Figure 16 which is analogous to the general 
diagram of Figure 1 and wherein: 

20 Ye(co) designates the input admittance of the resonator; 

MA is, in the case of an instrument where the string is 
struck, a model of hammer, expressing its speed from the 
strength f(t); 

Vs(t) is the speed of the string and Vh(t) that of the 
25 hammer; 

MV is a calculation model of the speed at the bridge which 
is then radiated by the sound board, from the strength and the 
speed of the string at the hammer-string contact point; 

G is the non-linear characteristic, and groups the non- 
30 linear function and the calculation means of the displacements 
Yh(t) and Ys(t) from Vh(t) and Vs(t); 

Vh(0) is the control parameter acting on the block MA, 
fixing the initial speed of the hammer at impact; 

L is the control parameter of the note played. 




In the case of an instrument where the string is rubbed, G is the 
non-linear friction characteristic, whereof numerous models can 
be found in the literature and whereof the control parameters 
are the pressure of the bow on the string and its speed of 
5 displacement. 

In a simplified model, the block MA may be deleted. 

In the case of an instrument where the string is struck, the 
discrete time model involves, as for certain elements of the 
acoustic models, the bilinear transform to approximate the time 
10 derivation operators. By noting W (noted V in the acoustic case) all 
the independent terms of (n) of the differential equation linking the 
speed of the string vs(n) and the strength f(n), and vh(n) the speed 
of the hammer, the transcription of the system of equations above 
in terms of sampled signals is: 

15 

f(n) = K(y h (n) - y B (n))*(l - c*(v h {n) - v a (n))) 

v h (n) = v h (n - 1) - _L_(/(n) + f(n - 1)) 

2J e Mh 

20 v.{n) = c Q f{n)+W (41) 
Whin) = y h (n - 1) + i(ifc(n) + v h (n - 1)) 

25 y s (n) = y s (n - 1) + 57-Mn) + v s (n - I)) 

Taking into account that the exponent (p) is not an 
integer, there is no explicit solution to this system, contrary to 
30 the acoustic case. By substituting in the first equation the 
expressions of yh(n) and ys(n) obtained from the two last 
equations, we obtain a « fixed point » type equation: 

35 /(») = (A - Bf(n)) p (C - Df(n)) 

Although there is no analytical solution, this type of 
equation is solved conventionally using iterative methods, the 
easiest being the so-called fixed point method. The regularity in 
the time domain of the strength f(n) enables to obtain very rapid 
40 convergence. Indeed, we observe that two or three iterations 
are sufficient. 




The model is controlled, for the note played, by acting on 
the resonator (length, diameter, tension of the string). The 
dynamics are controlled by the initial speed of the hammer, 
obtained by fixing vh(n=0). 

5 For exemplification purposes of realisation, Figure 17 

represents, in relation to time, at the top the speed of the string 
at the contact point (the eighth of its length), at the bottom the 
strength exerted by the hammer on the string, solutions of the 
previous system of equations, solved by the fixed point method. 

10 Similarly, Figure 18 represents the time trajectory of the 

strength in relation to the relative displacement of the hammer 
with respect to the string. 

We have thus described in detail an application of the 
method according to the invention to the simulation of a string 

15 instrument. As for the wind instruments described previously 
and contrary to the methods known previously, the invention 
enables to dispense with away-waves and return-wave 
quantities. 

Besides, it should be noted that the simulation model of a 
20 string instrument, illustrated in Figure 16 is very similar to the 
wind instrument model illustrated on Figure 1. Indeed, in both 
cases, they involve linear filters including delays, to realise 
non-linear interaction between two physical variables, so-called 
Kirchhoff variables, representative of the effect and of the 
25 cause of the phenomenon to be simulated. 

It appears therefore that the invention may be generally 
extrapolated to the simulation of any instrument operating by 
non-linear coupling between an excitation source and a 
resonator. 

30 For exemplification purposes, Figure 19 represents the 

general diagram of the model of such a digital instrument 
comprising, as usual, a control element I, a modelization 
element II and an element creating the sound III. 

As previously, the modelization element II includes a 

35 linear part 3 with a computing block (31) whose transfer function 




is, according to the instrument to be simulated, either the input 
impedance of the resonator Ze(co) or the admittance Ye(a>) and 
one non-linear part 2 which implements a non-linear function 
21. 

5 The block 1 may be a gestural sensor supplying control 

parameters CL acting on the linear part 3 of the model, and 
control parameters CNL acting on the non-linear part 2. 

According to the direction of the arrows indicated on the 
block 31 on Figure 19, the linear part 3 receives from the non- 

10 linear part 2, from left to right, when the transfer function of the 
computing block 31 is the impedance, an effect signal E to 
produce a cause signal C which is transmitted to the non-linear 
part 2, the latter producing, from this cause signal C, a new 
effect signal E bound for the linear part 3. 

15 Conversely, when the transfer function of the computing 

block 31 is the admittance, the linear part 3 receives from right 
to left, from the non-linear part 2, a cause signal C and 
produces an effect signal E which is transmitted to the non- 
linear part 2 to produce a new cause signal C bound for the 

20 linear part 3. 

The non-linear part 2 is associated with exciters 23 
transforming respectively the cause and effect signals to 
produce the other variables involved in the non-linear 
characteristic H. 

25 The block 4 includes calculation means of the sound to be 

emitted from the cause C and effect E signals, which is 
transmitted to a digital/analogue converter 5. 

The invention thus enables to simulate all sorts of 
instrument and is not limited, besides, to the field of music. 

30 Indeed, the method according to the invention could be 

also applied to the simulation of other oscillating phenomena, 
thanks to an adaptation of certain differential equations and a 
selection of other non-linear characteristics and of control 
parameters taking into account the physical characteristics of 

35 the phenomena to be simulated. 



